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Abstract 


The effect of switching between nonstiff and 
stiff methods on the efficiency of algorithms for 
Integrating chemical kinetic rate equations is 
presented. Different integration methods are 
tested by application of the packaged code LSODE 
to four practical combustion kinetics problems* 

The problems describe adl abati c t homogeneous gas- 
phase combustion reactions. It is shown that 
selective use of nonstiff and stiff methods in dif- 
ferent regimes of a typical batch combustion prob- 
lem Is faster than the use of either method for the 
entire problem. The Implications of this result to 
the development of fast integration techniques for 
combustion kinetic rate equations are discussed. 

Introduction 

The ordinary differential equations (ODE's) 
describing complex chemical reactions are char- 
acterized by widely different time constants. 
Although the differential equations are stable* 
standard numerical techniques such as the explicit 
Runge-Kutta and Adams methods are prohibitively 
expensive to use because of the severe steplength 
restriction Imposed by the requirements for numer- 
ical stability. Such systems of differential 
equations are commonly referred to as "stiff" 
systems, 

The problem of stiffness has been recognized 
for some time, e.g., 5 and several techniques have 
been developed for stiff ODE * s . ,At the present 
time, the packaged codes EPISODE 7 and LSODE 8 * 9 rep- 
resent the most extensively documented, tested and 
used routines for stiff ODE's. Amono several codes 
examined in recent detailed studies, 10-1 ^ LSODE was 
found to be the fastest for solving chemical kinet- 
ic rate equations, However, it Is recognized by 
combustion device modelers that LSODE is not fast 
enough for economical calculations of multidimen- 
sional reacting flow problems,* 5 

The numerical solution of combustion kinetic 
rate equations Is complicated by the existence of 
a narrow region ("heat release" zone) where the 
species concentrations and temperature change rap- 
idly, as illustrated in Fig. 1 for a typical batch 
reaction combustion problem. In the heat release 
regime, especially In the early part, many of the 
species and the temperature have positive time 
constants — an Indication that the governing ODE's 
are unstable. Since small steplengths are required 
for solving unstable ODE's, the use of methods 
designed for stiff problems — designated herein 
as "stiff methods" — may be inefficient. During 
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early heat release explicit "nonstiff methods" — 
i.e., methods suitable for nonstiff problems — 
may be adequate. However, Implicit methods are 
more accurate than explicit methods, which are 
therefore used only as predictors in predictor- 
corrector algorithms, 15 The corrector equations 
are Iterated until convergence Is obtained, It Is 
not clear what corrector Iteration technique is 
optimal in. the.nonstlff regime. Both simple or 
functional 16 * 17 and Jacobi-Newton 18 * 19 iteration 
techniques have been used because they avoid the 
expense associated with forming and inverting 
Jacobian matrices, which is required by Newton- 
Raphson iteration. However, much larger step- 
lengths can be used with Newton-Raphsbn iteration. 
For unstable ODE's this advantage may not be of 
much help and it is therefore not apparent which 
technique is the most efficient, 

During late heat release and equilibration the 
governing ODE's are stable so that Newton-Raphson 
iteration is the optimal convergence method. In 
these regimes, especially during equilibration, the 
different species approach the equilibrium state 
at different rates and the ODE * s are stiff — I.e., 
classical numerical techniques will require prohib- 
itive amounts of computer time In these regimes. 
Here, stiff methods are better suited to solving 
the problem, 

In developing an efficient algorithm to solve 
combustion kinetic rate equations, It is important 
to recognize and accommodate the widely different 
characteristics of the three regimes (Induction, 
heat release and equilibration) encountered In a 
typical combustion problem. Such a situation where 
the problem changes character, occurs in other 
areas and schemes have been proposed for automatic 
switching between stiff and nonstiff methods, 15 

The objective of the present Investigation is 
to examine the nature of the ODE'S arising in com- 
bustion chemistry. In particular, we examine the 
effect of switching between stiff and nonstiff 
methods on the computational work required to solve 
combustion kinetic rate equations. We also examine 
the use of different corrector Iteration techniques 
with nonstiff methods, 

Governing Ordinary Differential Equations 

The first order ODE 1 s describing the time rate 
of change of species i(i = 1,NS) can be written as 

dn. 

eft - 65 ^ 

(t = 0) «• given (1) 

T(t = 0) * given 

where ni is the mole number of species i; NS is 
the total number of distinct species in the qas 
mixture; T is the temperature; and f-j is net 
rate of formation of species 1 due to all forward 
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and reverse reactions In which species 1 partici- 
pates. A more detailed description of the govern- 
ing OOE’s Is given In Refs. 12 and 13. 

The Initial value problem Is to solve the sys- 
tem of Eqs. (1) for the chemical composition and 
temperature at the end of a prescribed time Inter- 
val, given the Initial conditions and the reaction 
mechanism. All problems considered In the present 
study Involve only adiabatic, homogeneous, gas 
phase chemical reactions. Tne problems are, how- 
ever, of two types — constant pressure and con- 
stant density. The following conservation 
equations serve as algebraic constraints on the 
species rate equations 

Constant pressure: 

NS 

y* V'i D ,j 0 “ constant (2a) 

1»1 


method devised by Nordsieck 22 to provide an Ini- 
tial estimate of the solution. To correct this 
estimate a range of iteration formulas Is Included 
In LSODE. The methods and corrector techniques 
attempted In this study are examined briefly; 
details are available In Refs. 7, 23 and 24. 

The ODE's presented In the previous section 
can be generalized as follows 



or using vector notation 

i s = ,f(jf) (4b) 

where 

y^ » I « 1,NS 


Constant density: 


£ n (h - RT) o U n constant 
id 1 1 o 


(2b) 


y NS+l “ T 
N « NS + 1 


(5) 


and an underscore represents a vectory quantity. 


where is the molal-speclflc enthalpy of 
species i ; R Is the universal gas constant: H 0 
and U 0 are the mass-specific enthalpy and inter- 
nal energy, respectively, of the Ideal gas mixture. 
Time differentiation of Eqs, (2a) and (2b) provide 
the following ODE's for the temperature: 


Constant pressure: 


dT 

BF' 


NS 

■ £ 

■w 


(3a) 


Constant density: 


NS 


- E f, h, - RT) 

dT 1=1 1 1 

cfT " "FIS' n Z TT 

£ n f (c pl “ R) 

1=1 


(3b) 


where c p j Is the constant-pressure molal specific 
heat of species 1. 

Methods and Iteration Techniques Examined 


The techniques used In this study are 
step-by-step methods. They compute approximations 
ZnU.yi, n!i *\M) to the exact solution 
ZltJ(s yj(t n );f = 1,N) at discrete points In time 
t n (n = 1,2,...). Assuming that solutions in_i, 
Yn_ 2 , • •• have been obtained at times tn_i, 
mi— 2 » ••• t * ie methods used In LSODE to advance the 
solution (» yj -) to time t n Involve linear 
multlstep formulas of the type 




.N 


( 6 ) 


where h n (* t n - t«_i) is the size of the step- 
length to be attempted; y^ J* f^(y k n ) t is the 
approximation to the exact*derivativ£ y^(t n ){= 
flCyk( t n3> ; and Ki» K^, aj, and are associ- 
ated with the formula selected to solve the problem 
over this time step. For the implicit Adams method 
of order q, K\ = 1, K 2 » q - 1 and Eq. (6) becomes 


y 1,n” y i,n-l +1 'n t n 6 J y i,n-j 1 


( 7 ) 


The objective of the present investigation was 
to examine the effect on the computational speed of 
using stiff and nonstiff methods in different 
regimes of a typical combustion kinetics problem. 

To accomplish this objective the packaged code 
LS0DE°> 9 was used because it contains both stiff 
and nonstiff methods and switching between the two 
methods is relatively straightforward. The methods 
included In this package are a variable-step, 
variable-order implicit Adams method, suitable for 
nonstiff problems, and a variable-step, variable- 
order backward differentiation formula (BDF), 20 
suitable for stiff problems. These methods are 
among the most efficient currently available for 
nonstiff and stiff problems, respectively. 17 * 21 
Both techniques employ a standard explicit predic- 
tor formula — a Taylor series expansion using the 


For a BDF of order q, « q, Kg * 0 and Eq. (6) 
becomes 

y 1,n = a j> y 1,n-J + h nVl,n 1 ’ M < fl > 

Equations (7) and (8) can be written in the general 
form 

y 1,n“ ¥ 1,n + l VVl,n 

-*1.ri + Wl( y k,n> 1 = 1 > N 

where n contains previously computed informa- 
tion. In*vector notation, Eq. (9) becomes 
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4“ VVtW (10) 

All of the different corrector Iteration tech- 
niques used In the present study to solve Eq. (10) 
can be generalized by the recursive relation 

(U) 

where I is the identity matrix, the matrix J 
depends on the Iteration method selected, (m) and 
(iti+1) denote the iteration numbers is the 

result obtained by the predictor step. 

The choice J « 0, called successive substitu- 
tion, 4 simple or functional Iteration 16 * 1 ' and Jacobi 
iteration, 1 results in 

d mH> ■ !, * wft”) <m 

Equation (12) is very simple to use but this method 
converges only linearly. 4 In addition, for suc- 
cessful convergence the steplength h n may be 
restricted to very small values, 1 

Newton-Raphson (NR) iteration, on the other 
hand, converges quadratically and can use much 
larger steplengths than functional Iteration. 4 * ib 
For this method J Is the Jacobian matrix and the 
elements Jfk(i,k « 1,N) are given by 

J ik .^ M-m (13) 


For this iteration technique much computational 
work can be required in forming the Jacobian matrix 
and In performing the linear algebra necessary to 
solve Eq. (11), To reduce this computational work 
the iteration matrix is not updated at every itera- 
tion. For additional savings it is updated only 
when the solution to Eq. (11) does not converge, 
Hence the Iteration matrix is only accurate enough 
for the iterations to converge and the same matrix 
may be used over several steps. 

The Jacobi-Newton (JN) Iteration technique 12 * 19 
can be obtained from the MR Iteration method by 
neglecting all off-diagonal elements of the 
Jacobian matrix. Hence for this technique 


Ik = of 


0; k * 1 


»v 


(14) 


This technique is os simple as functional iteration 
In the sense that no matrix inversion is Involved. 
Also it converges faster than functional Iteration 
— better than linear but not quite quadratic, iy 

In summary, functional and JN iteration tech- 
niques require much less work per step than NR 
iteration but have to use smaller steplengths and 
converge at slower rates. For stable problems 
where the Jacobian changes slowly NR iteration is 
clearly the optimal method. For unstable regimes, 
however, where rapidly changing solutions may 
require Sequent updating of the Jacobian for suc- 
cess* convergence, simple or JN Iteration may be 
more efficient, which may also be the case when 


very accurate numerical solutions are required. 
Because any change in the steplength alters the 
iteration matrix, It is not economical to consider 
small changes In the steplength with NR itera- 
tion. On the other hand, simple and JN Iteration 
techniques can take advantage of even modest in- 
creases In the ste, length. JN iteration requires 
a little more work per step than simple iteration 
but it converges faster* Also, it can use step- 
lengths at least as large as those used by simple 
iteration. 10 The optimal corrector technique 
therefore depends on the nature of the problem, the 
basic method used and the accuracy required of the 
numerical solution, 

In LSO0E both the basic method and the correc- 
tor iteration technique are selected via a method 
flag, MF. If NR Iteration is employed, either the 
user can provide analytical expressions for the 
elements of the Jacobian matrix or the code will 
estimate these elements by finite-difference 
approximations. For JN Iteration, however, this 
option is not available and the code uses 
internally-generated finite-difference approxima- 
tions for the diagonal elements of the Jacobian 
matrix. For all results obtained with NR Iteration 
analytical Jacoblans were used. The basic methods 
and iteration techniques employed in the present 
study are summarized in Table I, together with the 
relevant values for MF, 

Test Problems 


Four practical combustion kinetics problems 
were used in the present study. All four cases 
described adiabatic, homogeneous, gas-phase, tran- 
sient, batch combustion reactions. Test problem 1 
described the ignition and subsequent combustion 
of a mixture of 33 percent carbon monoxide and 
67 percent hydrogen with 100 percent theoretical 
air at an initial temperature of 1000 K and pres- 
sure of 10 atm. It involved 12 reactions among 11 
species. Test problem 2, involving 30 reactions 
among 15 species, described the Ignition and subse- 
quent combustion of a stoichiometric hydrogen-air 
mixture at 2 atm and 1500 K initial temperature, 
Both test cases 1 and 2 were at constant pressure 
and are discussed In more detail in Ref. 12. Test 
problem 3, taken from Burcat and Radhakr1$bnan,* b 
described the ignition and subsequent combustion of 
a stoichiometric propene-oxygen-argon mixture at 
an Initial temperature and pressure of 1700 K and 
4 atm, respectively, This constant density test 
case consisted of 113 reactions among 31 species, 
The reaction mechanism and rate constants were 
taken from Westbrook and Pitz. Zb Test case 4, 
taken from Bittker and Scull in, 2/ was a lean 
methane-air ignition and combustion problem at a 
constant pressure of 1 atm and initial temperature 
of 1645 K. This test problem involved 58 reactions 
among 24 species. 

Figure 1 presents the variation with time of 
the chemical species mole fractions and temperature 
for test problem 1, The variation of the tempera- 
ture with the reaction time for all four test cases 
is shown in Fig, 2. All four test problems were 
solved over a time period of 1 ms. This reaction 
period encompassed all three combustion regimes 
(induction, heat release and equilibration) for 
test problems 1-3. Test case 4, however, Included 
the first two regimes (induction and heat release) 
but only the beginning of equilibration. 
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Results 

In tills section we present the effects on the 
computational work of using stiff and nonstiff 
methods in different regimes of a typical combus- 
tion kinetics problem. All results were obtained 
on the NASA Lewis Research Center's IBM 370/3033 
computer using single-precision accuracy. 

As Illustrated In Fig. 1 and discussed in the 
section INTRODUCTION, a typical combustion kinetics 
problem consists of three distinctly different 
regimes; induction, heat release and equilibra- 
tion. During induction and early heat release 
when many of the ODE'S have positive time con- 
stants, small steprleogths must be used to insure 
solution accuracy. 1Z * 13 In these reoimes non- 
stiff methods may be more efficient. 15 During 
late heat release and equilibration when the ODE's 
are more stable, much larger steplengths can be 
used and NB iteration is the optimal convergence 
method. 5 * 15 * 19 In these later regimes, especially 
during equilibration, the ODE * s are stiff so that 
stiff methods are appropriate. 

To investigate if it is more efficient to use 
a nonstiff method during induction and early heat 
release, the variation of the computer time with 
the reaction time was examined for all values of 
the method flag, MF {» 10,11,13, and 21 — see 
Table I), used in this study. Pure relative error 
control is appropriate for the problems employed 
in this study, 1Z However, it could not be used 
because many of the mole numbers had zero Initial 
values. A mixed relative and absolute error con- 
trol was therefore used. Sufficiently small values 
for the local absolute error tolerances for the 
species were used to make the error control sub- 
stantially relative for mole fractions greater than 
0.1 ppm. For temperature pure relative error con- 
trol was used. To ensure that a comparison of com- 
putational work was made among comparably accurate 
methods, the same values for the absolute error 
tolerances were used with all methods and corrector 
iteration techniques. For clarity In presentation, 
methods corresponding to method flag MF = 10, 11, 
13, and 21 will hereafter be designated as methods 
10, 11, 13, and 21, respectively. 

Figures 3 and 4 present the variation of the 
computer (i.e., CPU) time (in seconds) with the 
reaction time for test problem 1 using values fgr 
the local relative error tolerance (EPS) of 10~ z 
and 10~ 5 , respectively. For method 10 (implicit 
Adams with functional iteration) and EPS a 10~ z the 
CPU time required up to the onset of heat release 
(reaction time -9 us, see Figs. 1 and 2) exceeded 
that required by method 21 (BDF with NR iteration 
using an analytical Jacobian) to solve the complete 
problem (Fig. 3), For EPS « KT 5 , however, the CPU 
times required during Induction and early heat 
release were about the same for both methods 
(Fig. 4). For methods 11 (implicit Adams with NR 
iteration using an analytical Jacobian) and 13 
(implicit Adams with JN iteration using internally 
generated approximations for the diagonal elements 
of the Jacobian matrix), the CPU times required 
during induction and early heat release compared 
favorably with, or were less than, those required 
by method 21. Note, however, that the CPU time 
required by method 21 for the complete problem was 
less than that required by all the nonstiff meth- 
ods, indicating that the problem was stiff. 


The results *lven In Figs. 3 and 4 show that 
JN and NR iteration techniques are more efficient 
than functional Iteration in the nonstiff regime. 
These results Indicate also that the use of a non- 
stiff method during Induction and early heat re- 
lease and a stiff method for the remainder of the 
problem would be more efficient than using either 
method for the complete problem. To examine the 
effects of such a switch the following procedure 
was used. The code was run up to reaction time 
t » t 5w i tch with a nonstiff method. After every 
step successfully executed by the routine, the 
value of the time reached by the integrator was 
checked to ensure that It did not exceed t$ w j tell* 
If the time exceeded t SW { tchp the method was 
switched to 21 and the problem was run to comple- 
tion with the stiff method. Upon completion of the 
problem, the CPU time required to solve the problem 
was calculated, In addition, the following per- 
formance parameters which give an Indication of the 
computational work required to solve the problem 
were noted; total number of steps required to 
solve the problem (NSTEP), total number of func- 
tional (i.e., derivative) evaluations (NFE) and 
total number of Jacobian evaluations (NJE). 

Different values for t sw u C h were attempted 
and the value resulting In the least CPU time to 
solve the problem was obtained by a trial-and-error 
process. Since the objective of the present inves- 
tigation was only to determine if switching methods 
resulted in efficiency increases and if so, to 
Identify the optimal Iteration technique to be used 
in the nonstiff regime, no attempt was made to in- 
corporate automatic method selection procedures. 

Table II presents the minimal CPU times ob- 
tained for test problem 1 using the two-stage solu- 
tion scheme outlined above and different Iteration 
techniques in the nonstiff regime. In this table 
^switch Is the reaction time (in us) up to which 
the program was run with the nonstiff method and 
the indicated iteration technique. For values of 
reaction time t > tc^tch the solution was 
obtained with the stiff method 21. Also given in 
Table II is the computational work required by 
method 21 to solve the complete problem. For 
method 10 and EPS > 10*^ the CPU times required up 
to the onset of heat release exceeded those re- 
quired by method 21 to solve the complete problem. 
Therefore no switching was attempted for these 
values or EPS and method 10. For EPS « 1(T 5 , how- 
ever, the combination of methods 10 and 21 was 
about 20 percent faster than method 21 for the com- 
plete problem (Table II). Note that fewer steps 
and functional evaluations were required by the 
s*.1ff method, indicating that the average step- 
length was smaller for method 10, However, the 
use of method 10 during induction and early heat 
release resulted in significantly fewer Jacobian 
evaluations* This was due to (a) not computing 
the Jacobian in the Initial regime and (b) fewer 
Jacobian evaluations being required In the second 
regime because of the use of smaller steplengths. 

Ihe combinations of methods 11 and 21 and of 
13 and 21 resulted in decreased CPU times (i.e., 
relative to method 21 for the complete problem) for 
most of the error tolerances (Table II). Also, in 
all cases the combination of nonstiff and stiff 
methods was faster than using the nonstiff method 
for the complete problem. Note that the time at 
which methods had to be switched generally 
Increased with decreasing EPS, i.e., Increasing 


accuracy requirement* This Implies that when EPS 
is decreased, accuracy requirements control the 
step size for a longer time, When accuracy re- 
quirements, and not numerical solution stability 
requirements, controLthe size of the step, the 
problem Is not stlff,*^ 1 Hence, the time over 
which It was more efficient to use a nonstiff 
method Increased with decreasing EPS. 

The combination of methods 11 and 21 resulted 
In CPU time decreases ranging from negligibly small 
to over 30 percent for test problem 1 (Table II). 
This switching process, l.e., use of a nonstiff 
method during Induction and early heat release and 
of a stiff method for the remainder of the problem, 
is not entirely satisfactory In that it does not 
always result In significant savings over the use 
of the stiff method 21 for the complete problem. 
Similar remarks apply to the use of method 13 In 
the Initial regimes. Note that for method 13 NJE 
includes two types of Jacobian matrix evaluations 
— the first number Is the total number of complete 
(l.e., analytical) Jacobian matrix evaluations re- 
quired and the second number Is the total number 
of diagonal matrix approximations (Table II). One 
difficulty encountered with the use of method 13 
was that it returned inaccurate solutions when 
relatively large values of EPS were used. This 
problem has been reported by others. It Is 
not clear If this was caused by poor approximations 
for the diagonal elements or by an unreliable con- 
vergence test. Another difficulty encountered with 
this method was serious numerical Instability for 
some test problems and values of EPS. Because of 
these problems with method 13 it was not attempted 
with the other three test cases. 

For the other three test problems and most of 
the error tolerances used, the runs with method 10 
required more CPU time until the onset of heat 
release than method 21 for tne complete problem, 
(e.g., Fig. 5). Hence, method 10 was also not 
attempted in the nonstiff regime for test problems 

2 to 4. 

Tables III, IV and V present the effects of 
switching between methods 11 and 21 for test prob- 
lems 2, 3 and 4, respectively. For purposes of 
efficiency comparison, the computational work re- 
quired by method 21 for the complete problem is 
also given In these tables. The results for test 
problem 2 (Table III) were very similar to those 
obtained for test problem 1. The use of the two- 
region scheme resulted In efficiency increases for 
most of the error tolerances and, as EPS was de- 
creased, the switching had to be performed at later 
times. 

For test problem 3, however, no significant 
efficiency increases could be obtained by using 
the nonstiff method 11 during Induction and early 
heat release and then switching to the stiff method 
21. But significant efficiency Increases could be 
obtained by switching before the onset of heat re- 
lease (Table IV). Note that for EPS * 10’* switch- 
ing from method 11 to method 21 at t = 0.03 ns 
(for this problem heat release started at about 

3 M s, Fig. 2) resulted In a CPU time decrease of 
over 40 percent. For test problem 3, unlike test 
problems 1 and 2, the temperature dropped by a sig- 
nificant amount (-21 K) during induction. This 
decrease in temperature was diagnosed by the code 
as an indication of stiffness, especially when low 
values were used for EPS. Note the sharp increase 


In CPU time incurred by the nonstiff methods during 
induction (Fig, 5), 

Test problem 4 was also quite different from 
test problems 1 and 2, Although the temperature 
drop during Induction was not significant (less 
than 1 K), this problem was characterized by a 
fairly long ignition delay period (Fig, 2). In 
addition, when the temperature started to Increase 
(at t - 20 us) It did so gradually and not as 
rapidly as for the other problems. For example, 
at t » 100 us the temperature had risen by only 
about 10 K. Unlike the other three test problems, 
test problem 4 included only the beginning of the 
equilibration regime, A nonstiff method was there- 
fore expected to be more efficient for most of the 
problem. However, the results given In Table V 
show that for Increased efficiency switching had to 
be performed during induction. This Indicates that 
for test problem 4 also It was more efficient to 
use a stiff method during Induction, as Illustrated 
In Fig. 6 for EPS » IQ" 4 . Note the larae Increase 
in CPU time for method 11 between t » I and 20 us, 
For method 21 the CPU time showed a large Increase 
between approximately 300 and 350 ps, corresponding 
to the rapid increase in the temperature between 
these times (Fig, 2), In this Interval method 11 
was more efficient (Fig, 6) because accuracy re- 
quirements control the step size. The effect of 
using a nonstiff method In this Interval was 
examined as follows for EPS * 10’ % The program 
was run with method 11 up to 2.5 us and between 
300 and 350 us. At all other times method 21 was 
used, This resulted In a total CPU time require- 
ment of 7.6 s — which was significantly faster 
than both the simple switch performed earlier 
(l.e., two-stage solution scheme) and method 21 
for the complete problem (Table V). 

The results discussed above Indicate that the 
induction regime is not necessarily nonstiff so 
that the use of a nonstiff method In this regime 
does not guarantee minimal computational work. In 
this regime the use of either a stiff method or a 
combination of nonstiff and stiff methods may re- 
quire the least computational work. To test this 
hypothesis the following procedure was adopted. 

The program was run with the stiff method 21 until 
the onset of heat release and also during late 
heat release and equilibration. During early heat 
release, however, a nonstiff method was used, 

Table VI presents the minimal CPU time 
obtained for test problem 1 using the three-region 
solution scheme discussed above — all iteration 
techniques were attempted during early heat 
release. In this table t sw \ and t sw ? are the 
times at which the methods were switched from non- 
stiff to stiff and from stiff to nonstiff, respec- 
tively. Note that as EPS was increased t S w A had 
to be decreased because heat release was predicted 
to start at an earlier time. As discussed previ- 
ously t sw ? had to be Increased with decreasing 
EPS. A comparison of Tables II and VI shows that 
for almost all iteration techniques and error tol- 
erances the three-stage solution scheme was faster 
than both the two-stage solution scheme proposed 
earlier and the stiff method 21 for the complete 
problem. Note that the use of this combination of 
stiff and nonstiff methods has resulted In about a c 
50 percent reduction in the CPU time for EPS * 10“* 
and method 13 during early heat release. Although 
the use of method 10 also resulted In efficiency 
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Increases, a very low value of EPS (10 - 5) was re- 
quired for significant reductions In the CPU time 
(Table VI), The use of such low values of EPS Is 
wasteful, especially for multidimensional reacting 
flow calculations. 30 This Indicates that either ON 
or NR Iteration should be used during early heat 
release. For small values of EPS JN Iteration 
(method 13) Is more efficient, Dut for large 
values of EPS NR Iteration (method 11) Is superior 
(Tabic VI). 

The results presented above Indicate that for 
efficient solution of combustion kinetic rate equa- 
tions, nonstiff methods should be used during early 
heat release, However, It Is not clear If JN or NR 
Iteration should be used In this regime, For large 
values of the local error tolerance JN Iteration 
resulted In significant errors, This could be due 
to the approximations for the Jacobian elements 
used In LSODE., .No.such problem was encountered 
with CREKIO 1 - 2 * which employs JN Iteration but 
with an analytical Jacobian. This suggests that 
JN Iteration with an analytical Jacobian should be 
attempted during early heat release. During late 
heat release and equilibration, however, a stiff 
method should be used. During induction either a 
stiff method or a combination of nonstiff and stiff 
methods appears to be the optimal choice. 

Conclusions 

A major conclusion of the present work Is that 
the combination of a nonstiff method during Induc- 
tion and early heat release and a stiff method dur- 
ing late heat release and equilibration does not 
always result In the optimal algorithm for solving 
combustion kinetic rate equations. During Induc- 
tion the use of either a stiff method or the combi- 
nation of nonstiff and stiff methods is Indicated. 
During early heat release a nonstiff method should 
be employed. However, It Is not evident If Newton- 
Raphson or Jacobl-Newton Iteration is the optimal 
convergence technique In the nonstiff regime. For 
large values of the local relative error tolerance 
the Jacobl-Newton iteration technique Included In 
the packaged code LSODE produced large errors and 
also resulted in unstable solutions. This may be 
the result of poor approximations for the Jacobian. 
Further experimentation, especially with an analyt- 
ical Jacobian, is necessary to resolve the question 
of which Iteration technique to select. During 
late heat release and equilibration stiff methods 
are optimal . 
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TABLE l, - SUMMARY OF METHODS AND CORRECTOR 
ITERATION TECHNIQUES EXAMINED 


Method 

flag* 

MF 

Basic method 

Iteration technique 

10 

11 

13 

Variable-step, 
variable order 
implicit Adams 

Simple or functional 
Nowton-Raphson with 
analytical Jacobian 
Jacobl-Newton with 
finite difference 
generated Jacobian 

21 

Variable-step, 
variable order, 
backward dif- 
ferentiation 
formula 

Mewton-Raphson with 
analytical Jacobian 


TABLE II. - SUMMARY OF COMPUTATIONAL WORK 
REQUIRED BY TWO-REGION SOLUTION FOR 
TEST PROBLEM 1 


Method 

EPS 

Switch* 

us 

NSTEP 

NFE 

NJE 

CPU, 

s 

10/21 

10-5 

17 

1701 

2801 

51 

3,87 

11/21 

io"| 

15 

99 

166 

30 

0,40 


IQ"? 

21 

173 

297 

33 

,63 


10-J 

23 

336 

540 

66 

1.17 


IO" 5 

40 

923 

1675 

139 

3,33 

13/21 

10-2 

13 

157 

308 

al6;45 

0.49 


io- 3 

14 

244 

459 

14; 60 

,64 


10-4 

13 

466 

958 

24 ; 135 

1.29 


IO' 5 

15 

956 

1742 

40; 146 

2,52 

21 

10-2 


115 

183 

30 

0,44 


IO" 3 

«... 

207 

346 

46 

.78 


10-4 



308 

504 

57 

1.08 


; IO" 5 

— 

1263 

2429 

255 

4.95 


a For Method 13 the first number Is the total 
number of complete Jacobian matrix evaluations 
and the second number Is the total number of 
diagonal matrix approximations. 





TABLE III. - SUMMARY OF COMPUTATION WORK 
REQUIRED BY TWO-REGION SOLUTION I*OR 
TEST PROBLEM 2 


Method 

EPS 

■Hi 

NSTEP 

m 

NJE 

m 

11/21 

10-2 

10-3 

ioi 

10- 5 

3 

5 

4.5 

20 

96 

159 

237 

2846 

158 

255 

368 

4686 

29 

43 

43 

422 

0.70 

1.07 

1.38 

16.5 

21 

10-2 

lO"? 

io*2 

icr 5 


100 

157 

295 

3527 

163 

243 

471 

5705 

29 

36 

63 

579 

0.71 

1.03 

1.87 

20.5 


TABLE IV. - SUMMARY OF COMPUTATIONAL WORK 
REQUIRED BY TWO-REGION SOLUTION FOR 
TEST PROBLEM 3 





































TABLE V. - SUMMARY OF COMPUTATIONAL WORK 
REQUIRED BY TWO-REGION SOLUTION FOR 
TEST PROBLEM A 


Method 

EPS 

^swl tch * 
us 

NSTEP 

HFE 

NOE 

CPU, 

5 

11/21 

icr\ 

25 

157 

254 

45 

2.24 


10-3 

12.5 

302 

521 

60 

4.14 


io*2 

2.5 

673 

1176 

192 

9 > 70 


lO" 5 

250 

1530 

2543 

210 

16.6 

21 

icr| 

r 1 

198 

326 

64 

2.93 


10-3 , 


490 

660 

172 

7.57 


l0 -4 1 


723 

1195 

207 

10.0 

i 

10“* , 

— 

2000 

i 

3708 

452 

27.1 


TABLE VI. - SUMMARY OF COMPUTATIONAL WORK REQUIRED BY 
THREE-REGION SOLUTION FOR TEST PROBLEM 1 



EPS 


tsw.2* 

MS 

NSTEP 

NFE 



NJE 

CPU, 

s 

21/10/21 

10-2 


9.0 

mm 

168 

31 

0.42 


10-3 

9.0 

10.0 

mm 

349 

37 

.71 


io-2 

9*5 

12.5 


633 

39 

1.09 


10* 5 

9.0 

14.0 

1076 

1726 

70 

2.81 



8,6 

15 

104 

157 

30 


10 -J 

8.5 

22 

167 

263 

34 

■Ell 


10-4 

9.0 

55 

278 

435 

43 

.90 


BM 

9.0 

65 

866 

1534 

132 

2.98 

21/13/21 


8.0 

13.0 

136 

231 

a 23;24 


10-3 

8.0 

13.5 

230 

403 

32;41 



10-4 

KIM 

13.5 

297 

499 

32 ; 33 



10-5 

ESS 

25 

981 

1691 

5B;153 



a For Method 13 the first nuiiiber is the total number of 
complete Jacobian matrix evaluations and the second 
number is the total number of diagonal matrix approxi- 
mations. 
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